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1  Introduction 


In  addition  to  requirements  such  as  strength,  stiffness,  and  aeroelastic  response  it  is  necessary 
to  design  both  metallic  and  composite  aircraft  structures  to  withstand  the  effects  of  damage. 
By  doing  so,  safe,  economical  fleets  with  high  operational  readiness  can  be  insured.  The 
importance  of  safety  and  operational  readiness  to  profitability  and  competitiveness  in  a 
commercial  aviation  setting  and  national  defense  in  a  military  one  is  apparent.  To  insure 
that  metallic  aircraft  structures  are  designed  and  maintained  to  withstand  the  effects  of 
damage,  the  Federal  Aviation  Administration  (FAA)  and  United  States  Air  Force  (USAF) 
have  established  specific  guidelines  which  must  be  followed.  The  FAA  requires  commercial 
transport  aircraft  certified  under  part  25  of  the  Federal  Aviation  Regulations  (FAR)  to  meet 
certain  Damage  Tolerance  Requirements  (DTR).  Similarly,  the  USAF  has  specified  a  set  of 
DTR  for  metallic  structures  which  are  described  in  detail  in  MIL-STD-1530A.  For  composite 
structures,  the  FAA  does  not  have  formal  requirements  for  certification  at  this  time.  The 
Air  Force,  on  the  other  hand,  has  specified  a  set  of  DTR  for  composite  structures  in  AFGS- 
87221A.  The  DTR,  as  set  forth  by  both  the  FAA  and  the  Air  Force,  essentially  state  that: 

•  the  residual  strength  of  an  airframe  structural  component  shall  not  drop  below  that 
required  to  sustain  limit  load,  and 

•  that  inspections  must  be  scheduled  to  insure  that  the  required  level  of  residual  strength 
is  maintained. 

The  primary  means  by  which  the  DTR  are  satisfied  by  commercial  and  military  air¬ 
frame  manufacturers  and  maintenance  organizations  is  through  the  performance  of  a  Damage 
Tolerance  Assessment  (DTA) .  The  DTA  involves  the  development  of  damage  growth  curves 
and  residual  strength  diagrams  for  individual  structural  components  of  an  airframe.  This 
allows  the  residual  strength  as  a  function  of  aircraft  usage  to  be  determined.  With  this  in¬ 
formation,  manufacturers  can  design  components  which  minimize  the  evolution  and  growth 
of  damage  and  its  effect  on  the  integrity  of  the  aircraft  structure.  Further,  appropriate  in¬ 
spection  intervals  can  be  specified  which  insure  that  the  structural  integrity  of  the  aircraft 
will  be  maintained  through  out  its  life. 

At  present,  there  is  no  software  which  integrates  damage  tolerance  with  the  other 
disciplines  (i.e.  strength,  stiffness,  aeroelastic  response,  etc.)  which  impact  the  design  of 
an  aircraft  structure.  As  a  result,  the  design  of  an  aircraft  structure  which  satisfies  damage 
tolerance  requirements  in  addition  to  those  of  other  disciplines,  is  presently  accomplished 
via  a  manual  “cut  and  try”  procedure.  This  type  of  design  process  is  time  consuming  and 
therefore  very  costly  [Nees  (1995)].  To  address  this  deficiency  in  existing  software,  this  SBIR 
project  is  implementing  a  damage  tolerance  module  into  the  multidisciplinary  analysis  and 
design  software,  ASTROS. 

Phase  I  of  this  SBIR  project  has  addressed  the  damage  tolerance  analysis  of  aircraft 
structures  made  of  laminated  composites.  The  damage  considered  was  in  the  form  of  a 
delamination  between  lamina  in  a  stiffened  panel.  The  Phase  I  effort  accomplished  the 
following. 
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•  A  damage  tolerance  software  applicable  to  laminated  composites,  called  BUCKDEL, 
was  developed.  The  BUCKDEL  software,  performs  geometrically  nonlinear  analysis 
of  stiffened  laminated  composite  plates  with  and  without  delaminations.  BUCKDEL 
allows  the  user  to  perform:  a  linear  static  solution;  a  linear  buckling  (eigenvalue)  anal¬ 
ysis;  and  a  nonlinear  post-buckling  analysis  through  both  limit  and  bifurcation  points. 
BUCKDEL  also  calculates  the  pointwise  energy  release  rates  around  a  delamination 
front  using  an  Equivalent  Domain  Integral  (EDI)  technique.  The  energy  release  rates 
can  be  used  to  predict  the  growth  and  onset  of  unstable  propagation  of  a  delamination. 

•  The  feasibility  of  using  a  global-local  approach  to  link  an  ASTROS  finite  element 
analysis  (global)  with  a  local  damage  analysis  (such  as  that  performed  by  BUCKDEL) 
was  investigated.  This  global-local  approach  was  found  to  be  workable,  due  in  large 
part  to  the  ASTROS  system  architecture  which  allows  the  user  to  introduce  special 
purpose  modules  by  making  use  of  the  SYSGEN  program  which  is  provided  with 
ASTROS. 

Based  on  the  Phase  I  results,  the  implementation  of  a  damage  tolerance  module  into  ASTROS 
which  accounts  for  typical  damage  in  both  composite  and  metallic  structure  is  feasible.  The 
remainder  of  this  report  will  describe  how  the  implementation  of  the  damage  tolerance 
module  will  be  carried  out.  In  addition.  Users  and  Theory  Manuals  for  BUCKDEL  are 
included  with  this  report.  The  BUCKDEL  software  represents  a  component  of  the  damage 
tolerance  module. 


2  SBIR  Technical  Objectives 


The  objective  of  this  SBIR  project  is  to  develop  damage  tolerance  software  which  can  be 
used  in  a  ‘stand  alone’  mode  or  as  an  analysis  module  in  the  multidisciplinary  analysis  and 
design  software,  ASTROS.  The  software  will  use  state  of  the  art,  computationally  efficient 
algorithms  for  determining  the  residual  strength  and  life  [Atluri  (1995)]  of  metallic  and  com¬ 
posite  structures  with  damage.  When  used  as  an  analysis  module,  it  will  enhance  the  existing 
capabilities  of  ASTROS  by  allowing  constraints  based  on  damage  tolerance  requirements  to 
be  considered  simultaneously  with  those  based  on  strength,  stiffness,  aeroelastic  response, 
etc.  during  design  optimization.  Included  in  potential  commercial  applications  of  such  a 
capability  are  industries  involved  in  the  design  of  aircraft  structures,  automobiles,  bridges 
and  buildings. 

ASTROS  is  well  suited  for  modeling  global  strength,  stiffness,  and  aeroelastic  response 
of  undamaged,  stiffened  structure.  However,  it  presently  does  not  have  the  capability  to 
account  for  local  damage  such  as  cracks,  delaminations,  and  penetration  holes.  Therefore, 
this  SBIR  project  will: 

•  utilize  ASTROS  existing  capabilities  for  performing  global  level  modeling  of  the  struc¬ 
ture; 
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•  utilize  an  existing  loads  program  called  USAGE  [Nees  (1995)]  to  general  load  spectra 
in  terms  of  ASTROS  load  cases; 

•  develop  a  finite  element  based  damage  tolerance  module  to  model  local  damage  in 
metallic  and  composite  materials; 

•  develop  interface  software  to  facilitate  the  linking  of  ASTROS  with  the  damage  toler¬ 
ance  module  via  the  SYSGEN  program  and  a  modified  MAPOL  sequence. 

•  define  new  bulk  data  entries,  relational  schema  and  error  messages  needed  to  integrate 
the  damage  tolerance  module  with  ASTROS. 

To  supplement  the  finite  element  based  damage  tolerance  module,  the  SBIR  project  will  also 
develop  ‘look  up’  tables  containing  cataloged  solutions  (i.e.  stress  intensity  factors,  energy 
release  rates,  etc.),  which  can  be  used  to  obtain  estimates  of  the  residual  strength  and  life 
of  common  structural  geometries.  In  addition  to  saving  computational  time,  these  ‘look  up’ 
tables  can  be  used  to  provide  sanity  checks  on  the  finite  element  calculations. 

The  damage  scenarios  that  will  be  considered  in  the  damage  tolerance  module  are: 

•  single  or  multiple  skin  cracks  (Widespread  Fatigue  Damage)  in  a  stiffened  structure, 
including  the  effect  of  broken  stiffeners  (Fig.  1)  [tensile  loading,  metallic  and  composite 
materials]; 

•  skin  crack  turning  at  stiffeners  and  its  effect  on  fail  safety  [tensile  loading,  metallic  and 
composite  materials]; 

•  single  or  multiple  part  elliptical  surface  flaws  at  holes  or  other  stress  raisers  (Fig.  2) 
[tensile  loading;  metallic  materials] 

•  arbitrary  shaped  holes  (Fig.  3) [compressive  loading,  metallic  and  composite  materials]; 
and 

•  arbitrary  shaped  delaminations  (Fig.  4) [compressive  loading,  composite  materials]. 

For  each  damage  scenario,  the  ability  to  calculate  the  residual  strength  of  a  structure  con¬ 
taining  that  type  of  damage  will  be  provided  in  the  damage  tolerance  module.  A  residual 
life  (fatigue)  capability  will  be  provided  for  elliptical  surface  flaws  and  skin  cracks  which 
propagate  in  a  self  similar  manner.  The  damage  tolerance  module  will  also  evaluate  DTR 
imposed  design  constraints  and  constraint  sensitivities  for  use  by  the  optimization  routine 
in  ASTROS. 

For  metallic  or  composite  aircraft  structure  loaded  in  tension,  the  damage  of  principal 
interest  during  design  is  usually  in  the  form  of  cracks  normal  to  the  direction  of  principal 
tension.  These  cracks  typically  occur  over  time  due  to  fatigue  or  suddenly  due  to  an  event 
such  as  an  uncontained  failure  of  a  rotating  engine  component  or  battle  damage.  A  crack 
arising  suddenly  due  to  a  catastrophic  event  is  an  example  of  discrete  source  damage  (DSD). 
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Figure  1:  Central  Crack  and  Broken  Stiffener  in  a  Panel 


Figure  2:  Failed  Lower  Wing  Panel  of  a  U.S.  Air  Force  C-141B  Due  to  the  Growth  of  a 
Surface  Flaw 


Figure  3:  Hole  in  the  Upper  (Compression)  Skin  of  a  Wing 


Figure  4:  Delamination  Due  to  Impact  of  a  Laminated  Composite 


In  reality,  DSD  in  a  structure  loaded  in  tension  such  as  a  fuselage  or  lower  wing  would 
be  in  the  form  of  an  irregular  shaped  hole.  However,  since  a  crack  of  length  ’a’  is  more 
critical  than  a  hole  of  diameter  ’o’,  the  crack  representation  of  DSD  represents  a  worst  case 
scenario.  Therefore,  for  reasons  of  conservatism  and  practicality  (both  experimental  and 
analytical)  the  crack  is  used  in  the  certification  of  primary  structural  elements  loaded  in 
tension.  Residual  strength  and  life  calculations  for  structure  containing  cracks  and  loaded 
in  tension  will  be  based  on  Linear  Elastic  Fracture  Mechanics  (LEFM).  Thus,  the  parts  of 
the  damage  tolerance  module  which  address  cracked  structure  will  compute  stress  intensity 
factors  and  their  sensitivities  to  changes  in  the  design  variables.  These  values  can  then 
be  used  to  evaluate  constraints  based  on  residual  strength  requirements.  For  constraints 
based  on  residual  life  (fatigue)  requirements,  computed  stress  intensity  factors  will  be  used 
in  crack  growth  equations  to  determine  the  time  required  for  a  crack  or  cracks  to  grow  to 
a  critical  size.  The  critical  size  is  determined  from  a  residual  strength  requirement  such  as 
that  requiring  the  structure  to  be  able  to  sustain  limit  load  at  any  time  in  its  life. 

For  metallic  or  composite  aircraft  structure  loaded  in  compression,  the  primary  concern 
is  with  DSD  in  the  form  of  a  hole  or  crack  parallel  to  the  direction  of  principal  compression. 
For  laminated  composites,  delaminations  between  lamina  are  also  of  concern.  Delaminations 
can  result  from  excessive  interlaminar  shear  stresses  or  through-the-thickness  tensile  stresses 
at  holes,  free  edges,  section  changes,  or  in  bonded  joints;  panel  buckling;  and  low  velocity 
impact.  The  modeling  of  delaminations  in  laminated  composites  was  addressed  in  Phase  I 
and  resulted  in  the  development  of  the  BUCKDEL  software.  The  residual  strength  of  struc¬ 
ture  loaded  in  compression  will  be  based  on  the  buckling  and  post-buckling  behavior  of  the 
structure.  For  laminated  composites  containing  delaminations,  the  pointwise  energy  release 
rate  around  the  delamination  front  will  also  be  used  in  the  residual  strength  prediction. 

The  motivation  for  modeling  DSD  in  a  structure  loaded  in  compression  in  the  form  of 
a  hole  or  crack  parallel  to  the  direction  of  principal  compression  is  as  follows.  For  a  stiffened 
structure,  the  buckling  load  is  expected  to  vary  significantly  with  the  size,  shape  (i.e.  circular 
or  elliptical),  and  location  (i.e.  distance  from  the  wing  tip)  of  the  hole.  Some  results  reported 
in  the  literature  [Vellaichamy  et.  al  (1990),  Nemeth  (1990),  and  Britt  (1994)]  indicate,  as 
is  to  be  expected,  that  for  an  elliptical  hole  with  the  major  axis  along  the  direction  of 
compression,  the  initial  buckling  load  is  lower  than  that  for  a  circular  hole  of  the  same  area. 
Similarly,  if  a  panel  is  subject  to  pure  shear  in  the  x-y  plane,  the  shear  buckling  load  will  be 
minimum  when  the  major  axis  of  the  ellipse  is  at  45  degrees  to  the  x  axis.  The  buckling  load 
will  decrease  with  an  increase  in  the  aspect  ratio  of  the  ellipse.  The  most  severe  case  being 
in  the  limit  when  the  ellipse  degenerates  into  a  crack  (with  the  crack  axis  along  the  direction 
parallel  to  the  direction  of  compression).  Based  upon  this  discussion,  it  is  anticipated  that 
the  worst  case  DSD  scenario  for  primary  structure  loaded  in  compression  will  be  a  crack 
oriented  such  that  its  axis  is  parallel  to  the  direction  of  maximum  compressive  stress. 
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Implementation  of  the  Damage  Tolerance  Module 


3.1  Damage  Tolerance  Analysis  of  Stiffened  Structure,  Using  a 
Global-Local  Methodology 

The  method  for  calculating  the  residual  strength  and  life  of  a  stiffened  structure  made  of 
metallic  and/or  composite  materials  is  to  use  a  global-local  methodology.  The  global  analysis 
will  be  accomplished  via  ASTROS  while  the  damage  tolerance  module  will  be  used  for  the 
local  analysis.  To  illustrate  this  methodology,  consider  the  case  of  a  wing  containing  a  crack 
in  the  skin  of  the  lower  surface  [Fig.  5].  (Note  that  the  methodology  is  general  in  nature 
and  can  be  used  for  structure  containing  damage  in  the  form  of  penetration  holes  and/or 
delaminations.)  A  coarse  finite  element  mesh  is  first  used  to  model  the  global  behavior  of 
the  cracked  wing.  The  traction  and/or  displacement  boundary  conditions  to  be  applied  to 
the  local  model  are  determined  from  the  results  of  the  global  analysis.  These  boundary 
conditions  include  the  reaction  forces,  exerted  by  the  stringers  and  ribs,  on  the  skin.  The 
Finite  Element  Alternating  Method  (FEAM)  is  used  in  the  local  analysis  to  determine  the 
stress  intensity  factors.  The  FEAM  allows  a  coarse  finite  element  mesh  to  be  used  for  the 
local  model  of  the  skin  because  the  crack  tip  fields  are  captured  by  an  analytical  solution  and 
thus  cracks  need  not  be  modeled  explicitly  in  the  mesh.  The  FEAM  along  with  the  other 
local  damage  modeling  methodologies  to  be  implemented  in  the  damage  tolerance  module 
will  be  discussed  in  detail  in  later  sections. 

In  the  global  analysis,  the  stringers  are  modeled  as  beams  (ASTROS  GEAR  elements) 
attached  to  the  skin  (ASTROS  CQUAD4  or  CTRIA3  elements),  and  ribs  are  modeled  as 
plates  (ASTROS  CQUAD4  or  CTRIA3  elements)  or  shear  panels  (ASTROS  CSHEAR  ele¬ 
ments),  as  illustrated  in  Fig.  6.  To  account  for  fastener  flexibility,  springs  (ASTROS  CELASl 
or  CELAS2  element)  can  be  used  to  model  connections  between  stiffening  elements  and  the 
skin.  This  level  of  detail  in  the  global  model  should  be  sufficient  for  most  situations.  How¬ 
ever,  the  user  is  free  to  either  decrease  or  increase  the  complexity  of  the  global  model,  within 
the  limits  imposed  by  the  modeling  capabilities  of  ASTROS. 

Damage  in  the  form  of  delaminations  typically  need  not  be  accounted  for  in  the  global 
model.  On  the  other  hand,  damage  in  the  form  of  cracks  or  penetration  holes  does  need  to  be 
accounted  for  if  it  is  of  suiSScient  size  to  affect  load  transfer  in  the  structure.  To  accomplish 
this,  holes  are  modeled  explicitly  (as  holes)  while  cracks  are  modeled  via  unconnected  nodes 
at  the  crack  locations.  This  crude  representation  of  the  crack  in  the  global  model  reflects  the 
loss  of  stiffness  in  the  structure,  so  that  the  redistribution  of  loads  among  the  skin,  stringers 
and  ribs  can  be  captured.  Broken  stringers  and  ribs,  if  any,  are  also  accounted  for  in  a  similar 
manner.  The  details  of  the  crack  tip  fields  are  ignored  at  this  level  of  analysis.  The  global 
analysis  results  are  used  to  construct  a  free-body  diagram  (Fig.  7)  of  the  cracked  sheet  (local 
model),  with  the  applied  loading  on  the  sheet  being  the  reaction  forces  from  stringers  and 
ribs  as  well  as  the  loading  on  the  periphery  of  the  sheet. 

For  problems  involving  fatigue  where  the  crack  is  of  sufficient  size  that  it  must  be 
accounted  for  in  the  globlal  model,  it  will  be  necessary  to  update  the  crack  size  periodically. 
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Figure  5:  Flow  Chart  Illustrating  the  Damage  Tolerance  Analysis  of  a  Wing  Using  a  Global- 
Local  Methodology 
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This  will  allow  changes  in  the  load  distribution  to  be  determined  and  the  boundary  conditions 
applied  to  the  local  model  to  be  updated.  In  general,  the  same  global  model  can  be  used  for 
a  number  of  local  analyses.  The  SBIR  project  will  automate  this  updating  of  crack  size  in 
the  global  model. 

During  the  course  of  the  Phase  I  project,  the  possibility  of  using  an  intermediate  analysis 
between  the  global  and  local  analysis  was  considered.  The  intermediate  analysis  would 
provide  an  additional  level  of  detail  which  in  some  cases,  may  improve  the  fidelity  of  the 
local  boundary  conditions.  The  global-local  approach  was  settled  upon  over  the  global- 
intermediate-local  approach  for  two  reasons.  First,  since  ASTROS  is  primarily  a  preliminary 
design  tool,  the  geometric  details  required  to  do  the  finite  element  modeling  of  stiffeners  and 
fasteners  at  the  intermediate  level  may  not  be  available  at  this  design  phase.  Second,  the 
additional  computational  burden  imposed  by  the  intermediate  analysis  would  likely  slow 
down  the  ASTROS  iterative  design  process.  This  would  be  of  particular  concern  during 
the  design  of  a  structure  such  as  a  wing  in  which  damage  related  constraints  were  being 
considered  in  multiple  panels. 

3.2  Damage  Tolerance  Module 

The  finite  element  based  damage  tolerance  module  for  modeling  local  damage  will  implement 
the  following  computational  methods: 

•  Finite  Element  Alternating  Method  [Wang  and  Atluri  (1996)]  for  modeling  single  or 
multiple  cracks  (including  Widespread  Fatigue  Damage)  under  tension  loading; 

•  hybrid  crack-tip  finite  elements  [Atluri  (1986)]  for  modeling  single  or  multiple  cracks 
under  tension  loading  and  crack  turning  at  a  stiffener; 

•  Finite  Element  Alternating  Method  for  modeling  single  or  multiple  part  elliptical  sur¬ 
face  flaws  at  holes  or  other  stress  raisers  [Nishioka  and  Atluri  (1983)];  and 

•  a  multi-domain  modeling  method  (implemented  in  BUCKDEL  during  Phase  I)  for 
calculating  the  buckling  and  post-buckling  strengths  of  stiffened  panels  with  arbitrary 
shaped  holes  and/or  delaminations. 

This  module,  when  combined  with  ASTROS  will  provide  a  global-local  methodology  for 
modeling  damage  in  a  complex  stiffened  structure.  This  will  allow  residual  strength  and 
residual  life  (fatigue)  constraints  stemming  from  the  DTR  to  be  treated  by  ASTROS  dur¬ 
ing  multidisciplinary  design  and  optimization.  In  addition  to  finite  element  based  damage 
modeling,  ‘look  up’  tables  containing  cataloged  solutions  will  provide  a  means  of  treating 
problems  involving  common  geometries  as  efficiently  as  possible  as  well  as  a  convenient  way 
of  checking  the  finite  element  calculations. 
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Figure  8:  Superposition  Principle  Used  in  the  Finite  Element  Alternating  Method 

3.2.1  Finite  Element  Alternating  Method  For  Modeling  Single  or  Multiple 
Cracks:  Tensile  Loading 

Once  the  global  analysis  is  performed  with  ASTROS,  one  can  consider  the  free-body  diagram 
of  the  cracked  skin  alone  (see  Fig.  7);  the  skin  being  subjected  to  far-field  tractions,  and  the 
stiffener  reaction  forces.  The  stress-intensity  factors  for  single  or  multiple  cracks  (including 
Widespread  Fatigue  Damage)  in  the  skin  can  be  determined  in  the  local  analysis  using  the 
Finite  Element  Alternating  Method  (FEAM),  while  still  using  a  coarse  finite  element  mesh. 
The  problem  in  Fig.  7  can  be  solved  with  the  FEAM  depicted  in  Fig.  8,  wherein  it  can  be 
seen  that  the  problem  of  Fig.  7,  can  be  identified  with  the  problem  labeled  as  the  “original 
problem”  in  Fig.  8. 

The  FEAM  is  based  on  the  Schwartz-Neumann  alternating  method  (see  the  Appendix 
for  theoretical  details).  The  combination  of  the  powerful  finite  element  method,  and  the 
analytical  solutions  for  cracks  in  an  infinite  domain,  subjected  to  arbitrary  crack  surface 
tractions  (that  have  only  recently  been  developed),  enable  the  alternating  method  to  be  used 
to  solve  fracture  problems  for  complex  structures.  With  the  help  of  the  initial  stress  method, 
it  can  even  be  used  in  elastoplastic  analyses  [Wang,  Brust  and  Atluri  (1995a),  (1995b),  and 
(1995c)].  It  has  been  used  successfully  in  the  evaluation  of  residual  strength  and  life  of 
cracked  metallic  and  composite  structures.  The  Finite  Element  Alternating  Method  solves 
the  problem  of  cracks  in  finite  bodies  by  iterating  between  the  analytical  solution  for  an 
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embedded  crack  in  an  infinite  domain,  and  the  finite  element  solution  for  the  uncracked 
finite  skin.  Essentially,  it  is  a  fixed  point  iteration  scheme  which  solves  the  superposition  of 
the  following  two  problems  [see  Fig.  8]. 

1.  the  uncracked,  finite-sized  skin  subjected  to  external  loads  (including  the  reaction 

forces  exerted  by  the  stiffeners  on  the  skin)  and  unknown  external  boundary  loads; 

2.  a  crack  in  an  infinite  sheet  subjected  to  a  crack  surface  traction 

The  crack  surface  traction  in  the  infinite  sheet  cancels  out  the  cohesive  traction  at  the 
location  of  the  crack  in  the  first  problem;  while  the  unknown  external  boundary  loads  in  the 
first  problem  cancel  out  the  tractions  at  the  same  location  in  the  second  problem.  Thus, 
the  original  problem,  i.e.  a  cracked  sheet  of  a  finite  dimension  subjected  to  reaction  forces 
exerted  by  the  stiffeners  with  other  boundary  loadings,  is  exactly  the  superposition  of  these 
two  problems. 

The  FEAM  is  used  to  solve  this  superposition  problem.  First,  reverse  the  cohesive 
tractions  at  the  location  of  the  crack  in  the  finite  element  model  of  the  uncracked  skin  and 
use  them  as  the  load  acting  on  the  crack  in  an  infinite  sheet.  Then,  reverse  the  residuals  at 
the  locations  of  the  far  field  boundries  in  the  infinite  sheet  and  apply  them  as  loads  acting  on 
the  boundaries  of  the  uncracked  skin.  In  this  way,  the  cohesive  tractions  at  the  location  of 
cracks  and  residuals  at  the  locations  of  the  external  boundaries  are  corrected  iteratively.  This 
procedure  converges  very  fast,  usually  in  two  or  three  iterations.  A  flow  chart  illustrating 
the  FEAM  is  shown  in  Fig.  9. 

Fracture  mechanics  parameters  can  be  found  accurately  because  the  near  crack  tip  fields 
are  captured  exactly  by  the  analytical  solutions.  Coarser  meshes  can  be  used  in  the  finite 
element  analysis  because  the  cracks  are  not  modeled  explicitly.  The  finite  element  method 
is  only  used  to  compute  the  cohesive  tractions  at  the  crack  location,  which  has  a  smooth 
distribution.  Therefore,  a  very  coarse  mesh  can  be  used.  Fig.  10  shows  the  typical  finite 
element  meshes  around  the  crack  tip,  when  a)  the  Equivalent  Domain  Integral  (EDI)  based 
method  is  used  to  evaluate  stress  intensity  factors;  or,  b)  the  FEAM  is  used.  In  Fig.  10,  the 
EDI  based  method  also  uses  singular  quarter-point  elements.  The  simplicity  of  the  FEAM 
mesh  relative  to  the  EDI  mesh,  which  must  explicitly  model  crack  tips  is  apparent  from  this 
figure. 

In  a  parametric  analysis  of  various  crack  sizes,  such  as  is  necessary  in  fatigue  calcu¬ 
lations,  the  stiffness  matrix  of  the  finite  element  model  is  decomposed  only  once,  since  the 
stiffness  of  the  uncracked  structure  remains  the  same  for  all  crack  sizes.  In  the  other  ap¬ 
proaches,  such  as  those  using  singular/hybrid  type  special  crack-tip  finite  elements  or  EDI 
methods,  the  cracks  must  be  modeled  explicitly.  Therefore,  the  global  stiffness  matrix  must 
be  computed  and  decomposed  for  each  crack  size.  Thus,  the  FEAM  is  very  efficient  in 
terms  of  both  computational  time  and  human  effort  (i.e.  mesh  generation)  when  applied  to 
problems  such  as  fatigue  crack  growth. 

Finally,  it  is  noted  that  a  simple  superposition  method  can  be  used  to  construct  the  so¬ 
lution  for  multiple  cracks  in  an  infinite  domain,  subjected  to  arbitrary  crack  surface  tractions. 
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Figure  9:  Flow  Chart  of  the  Finite  Element  Alternating  Method 
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Figure  10:  The  Finite  Element  Mesh  When  a)  The  EDI  Method  is  Used;  b)  The  Finite 
Element  Alternating  Method  is  Used 

using  the  solution  for  a  single  crack  in  an  infinite  domain  (see  the  Appendix  for  theoretical 
details).  With  the  solution  for  multiple  cracks  in  an  infinite  domain,  the  FEAM  can  be  used 
to  solve  problems  of  multiple  cracks  with  arbitrary  crack  lengths  and  orientations  at  arbi¬ 
trary  locations.  This  is  particularly  useful  in  the  treatment  of  Widespread  Fatigue  Damage 
(WFD). 

3.2.2  Hybrid  Crack-Tip  Finite  Elements  for  Modeling  Skin-Crack  Turning  at 
Stiffeners:  Tensile  Loading 

In  assessing  the  integrity  of  a  structure  containing  a  crack,  accurate  evaluation  of  fracture 
parameters  (i.e.  stress  intensity  factors)  is  required.  If  the  finite  element  method  is  used 
for  such  purposes,  proper  numerical  modeling  of  the  crack-tip  singularities  is  necessary.  An 
alternative  to  the  Finite  Element  Alternating  Method  is  to  use  singular/hybrid  crack-tip 
finite  elements.  The  use  of  these  special  elements  enables  one  to  use  a  relatively  coarse 
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finite  element  mesh  around  the  crack  tips;  as  compared  to  the  very  fine  and  focused  mesh 
when  ordinary,  non-singular  finite  elements  are  used.  However,  explicit  crack-tip  meshing 
must  still  be  carried  out  and  hence  this  method  requires  considerable  human  effort  in  the 
construction  of  meshes  as  compared  to  the  FEAM. 

The  hybrid  crack  element  can  be  used  to  analyze  general  anisotropic  materials.  Hybrid 
multilayer  elements  [Nishioka  and  Atluri(1980)]  can  be  used  to  study  complex,  non-self¬ 
similar  crack  growth  in  composite  structures.  While  self-similar  crack  growth  in  a  homo¬ 
geneous  anisotropic  sheet,  can  be  studied  easily  using  the  FEAM,  the  FEAM  cannot  be 
applied  to  non-self-similar  crack  growth.  Thus  the  hybrid  crack  element  will  be  used  in  the 
SBIR  project  to  handle  problems  involving  non-self-similar  crack  growth.  In  hybrid  multi¬ 
layer  elements,  the  fully  three  dimensional  stress-state,  including  <733  is  accounted  for.  The 
mixed-mode  stress  and  strain  singularities  whose  intensities  vary  within  each  layer  (or  a 
group  of  layers  in  a  repeated  lay-up  sequence)  near  the  crack-front,  are  built  into  the  formu¬ 
lation  a  priori.  The  interlayer  traction  reciprocity  conditions  are  satisfied  through  the  use 
of  Lagrange  multipliers  and  individual  cross-sectional  rotations  of  each  layer  are  allowed. 

In  aircraft  design,  the  problem  of  crack-turning  near  a  stiffener  is  of  particular  interest. 
This  is  desirable  in  most  situations  as  it  acts  as  a  form  of  crack  arrest.  In  fact,  most 
aircraft  structures  are  designed  with  this  fail  safe  feature.  The  analysis  of  crack  turning  at 
a  stiffener  will  be  handled  in  the  damage  tolerance  module  by  using  the  hybrid  crack-tip 
element  method.  If  a  crack  approaches  the  stiffener  at  right  angles,  the  T-stress  at  the  crack 
tip  plays  an  important  role  in  postulating  a  criterion  which  will  predict  whether  the  crack 
will  turn  along  the  stiffener.  During  the  process  of  crack-turning,  the  hybrid  crack  elements 
will  be  used  near  the  tip  of  the  curving  crack,  to  evaluate  the  condition  for  the  continued 
turning  of  the  crack.  A  constraint  based  on  the  T-stress  at  the  crack  tip  will  be  introduced 
in  order  to  account  for  this  crack-turning  criterion  in  the  ASTROS  design  iteration. 

3.2.3  Finite  Element  Alternating  Method  for  Modeling  Single  or  Multiple  Part 
Elliptical  Surface  Flaws  at  Holes  or  Other  Stress  Raisers:  Tensile  Loading 

The  Finite  Element  Alternating  Method  for  through  skin  cracks  has  already  been  described. 
The  application  of  the  FEAM  to  surface  cracks  is  similar.  In  the  FEAM  for  a  surface  crack 
in  a  finite  solid,  two  solutions  are  required.  Solution  1:  A  general  analytical  solution  for  an 
embedded  elliptical  crack  in  a  body  subject  to  arbitrary  crack  face  tractions  [Vijayakumar 
and  Atluri  (1981)]. 

Solution  2:  A  numerical  scheme  (the  finite  element  method  in  this  instance)  to  solve  for  the 
stresses  in  an  uncracked  finite  body.  The  analytical  solution  is  for  a  crack  of  elliptical  shape. 
Because  of  this,  all  physical  cracks  that  are  being  analyzed  must  be  either  elliptical  or  part 
elliptical  in  shape. 

The  finite  element  solution  involves  the  analysis  of  the  uncracked  solid.  Thus,  non-zero 
stresses  are  calculated  at  the  location  of  the  actual  crack.  These  stresses  must  be  removed 
in  order  to  create  a  traction  free  crack  as  in  the  actual  problem.  The  infinite  body  with  an 
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embedded  crack  has  a  solution  which  is  valid  for  an  arbitrary  distribution  of  tractions  on 
the  crack  face.  The  detailed  steps  involved  in  the  FEAM  for  a  crack  in  a  finite  body  are  as 
follows. 

1.  Solve  the  uncracked  finite  body  under  the  given  external  loads  using  the  finite  element 
method.  The  uncracked  body  has  the  same  geometry  as  in  the  given  problem  except 
for  the  crack.  For  example,  when  a  crack  emanates  from  a  hole  in  a  structure,  the  hole 
must  still  be  analyzed  in  the  uncracked  structure. 

2.  Using  the  finite  element  solution,  the  program  computes  the  stresses  at  the  crack 
location. 

3.  It  then  compares  the  residual  stresses  calculated  in  Step  2  with  a  permissible  stress 
magnitude. 

4.  The  residual  stresses  at  the  crack  location  as  computed  in  Step  2  are  reversed  to  create 
the  traction  free  crack  faces  as  in  the  given  problem.  Prom  this,  the  program  determines 
a  polynomial  form  for  these  stresses  using  a  ’’least  squares  fit”. 

5.  The  anal3d;ical  solution  to  the  infinite  body  problem  with  the  crack  subject  to  the 
polynomial  loading  calculated  in  Step  4  is  now  obtained. 

6.  The  stress  intensity  factors  for  the  current  iteration  are  then  calculated  from  the  ana¬ 
lytical  solution. 

7.  The  residual  stresses  on  the  external  surfaces  of  the  body  due  to  the  applied  loads  on 
the  crack  faces,  are  now  computed.  To  satisfy  the  given  traction  boundary  conditions 
at  the  external  boundaries,  the  residual  stresses  on  the  external  surfaces  of  the  body 
are  reversed  and  this  allows  calculation  of  the  equivalent  nodal  forces. 

8.  Consider  the  nodal  forces  in  Step  7  as  externally  applied  loads  acting  on  the  uncracked 
body. 

All  the  steps  in  the  iteration  process  are  repeated  until  the  residual  stresses  on  the  crack 
surface  become  negligible.  It  has  been  observed  that  this  iteration  process  typically  takes 
three  or  four  steps.  The  overall  stress  intensity  factor  solution  is  obtained  by  adding  the 
stress  intensity  factor  solutions  for  all  iterations. 

One  recent  application  [O’Donoghue,  Atluri  and  Pipkins  (1995)]  of  the  FEAM  was 
to  the  analysis  of  fatigue  cracking  in  the  lower  wing  skin  of  the  U.S.  Air  Force  C-141B 
(Fig.  11).  In  this  application,  the  growth  of  corner  cracks  from  the  weep  holes  located  in 
the  integral  risers  of  the  wing  skin  was  modeled  with  the  FEAM  (Fig.  12  and  Fig.  13). 
Fatigue  crack  growth  predictions  were  made  wherein  the  FEAM  was  used  to  generate  stress 
intensity  factors  which  were  inturn  used  in  a  Forman  equation.  The  load  spectrum  used  was 
comprised  of  peak-valley  pairs  representing  3027  equivalent  flight  hours.  Comparisons  with 
limited  test  data  showed  good  correlation  between  the  FEAM  based  fatigue  crack  growth 
predictions  and  the  experimental  data. 


16 


Figure  11:  U.S.  Air  Force  C-141B 


Figure  13:  Cross-Section  of  Failed  Lower  Wing  Panel  of  the  C-141B 

3.2.4  Multi-Domain  Modeling  of  the  Buckling  and  Post-Buckling  Strengths  of 
Stiffened  Panels  with  Arbitrary  Shaped  Holes  and  Delaminations:  Com¬ 
pressive  Loading 

The  primary  damage  in  structure  loaded  in  compression  which  must  be  considered  during 
design  is  a  hole  (Fig.  14).  For  laminated  composites,  delaminations  will  also  be  of  relevance. 
It  is  to  be  expected  [with  some  preliminary  results  that  confirm  these  expectations  being 
reported  in  Vellaichamy,  et.  al  (1990)]  that  the  buckling  load  of  a  panel,  with  a  hole,  in 
compression  will  depend  on  the  shape  and  the  size  of  the  hole.  The  buckling  load  of  a  panel 
with  an  elliptical  hole  whose  major  axis  is  aligned  with  the  direction  of  compression  can  be 
expected  to  be  lower  than  that  with  a  circular  hole  of  the  same  area.  Further,  it  is  expected 
that  for  an  elliptical  hole  of  a  given  area,  the  buckling  load  will  decrease  with  an  increase 
in  the  aspect  ratio  of  the  ellipse,  as  long  as  the  major  axis  of  the  ellipse  is  aligned  with  the 
direction  of  compression.  In  the  limit  as  the  elliptical  hole  shrinks  to  a  crack  whose  axis 
is  parallel  to  the  direction  of  compression,  with  the  slightest  imperfection,  the  crack  may 
propagate  in  mode  III  in  post-buckling  deformation.  Such  mode  III  behavior  would  severely 
affect  the  structural  integrity.  BUCKDEL,  as  developed  in  Phase  I  of  this  SBIR  project,  can 
be  used  to  compute  the  buckling  and  post-buckling  response  of  stiffened  structure  containing 
damage  in  the  form  of  holes  and  delaminations  of  arbitrary  shape.  In  addition,  it  computes 
the  pointwise  energy  release  rate  around  the  delamination  front.  In  future  work,  the  ability 
to  treat  the  mode  III  crack  problem  will  be  added  to  BUCKDEL. 

A  brief  overview  of  BUCKDEL  is  given  here.  For  more  details,  see  the  Users  and  Theory 
Manuals  which  are  included  with  this  report.  BUCKDEL  uses  a  multi-domain  method  to 
model  delaminations  of  arbitrary  shape.  In  this  method,  the  delaminated  shell  is  assumed 
to  be  assembled  with  three  regions-(l)  Undelaminate:  undelaminated  zone;  (2)  Delaminate: 
thinner  side  of  the  delaminated  zone  and  (3)  Base:  thicker  side  of  the  delaminated  zone. 
Transverse  shear  deformation  plays  an  important  part  in  the  case  of  composite  structures. 
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Figure  14:  Damage  in  the  Upper  (Compression)  Skin  of  a  Wing 


hence,  it  is  introduced  explicitly  and  the  assumptions  of  Reissner-Mindlin  theories  of  plate 
bending  are  used  for  modeling  each  region  of  the  multi-domain  model.  Thus,  for  each 
region,  the  3-dimensional  displacement  field  (U  =  {Ui  U2  U3})  is  expressed  in  terms  of  the 
corresponding  mid-surface  displacement  (u  s  {m  U2  U3})  and  rotation  {0  =  {^i  $2  0})  fields 
as, 

U«(^a, X3)  =  u«(a;,)  -  (1) 

where  (a  =  1, 2)  are  the  inplane  curvilinear  shell  coordinates  and  is  the  thickness 
coordinate  for  the  (i  =  1, 2, 3)  shell.  The  structural  continuity  at  the  delamination  front 
r  is  maintained  by  assuming  the  deformation  to  be  unique  at  the  junction  of  the  three  shells 
i.e.  on  F  in  accordance  with  the  Reissner-Mindlin  law  of  flexure  (Eq. 

(1)).  In  other  words,  at  the  delamination  edge,  the  mid-surface  degrees  of  freedom  of  the 
delaminate  and  the  base  shells  are  assumed  to  be  related  to  those  of  the  undelaminated  shell 
by, 
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where  is  the  distance  of  the  midsurface  of  the  i^^  shell  from  the  laminate  midsurface. 

Each  lamina  is  assumed  to  be  orthotropic  and  the  inplane  stresses  =  {an  022 
and  the  transverse  shear  stresses  =  {tiz  r23}(*)  are  related  to  linear  components  of 
membrane  strain  £22  £^12}^*^;  nonlinear  components  of  membrane  strain  = 

[uii  U22  ^’12}^*^;  flexural  strain  due  to  mid-surface  rotation  =  {/cn  K22  «i2}^*b  flexural 
strain  due  to  transverse  shear  strain  =  {xii  X22  Xi2}^*^  and  transverse  shear  strains 
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where  the  material  constitutive  terms  E^f  are  functions  of  the  thickness  coordinate  of  each 


shell  ^3^  Generally,  for  a  laminate  with  orthotropic  layers,  Ejf  are  assumed  to  be  piecewise 
constants  over  the  laminate  thickness. 

Integrating  along  the  thickness,  the  constitutive  equations  can  be  written  in  terms  of 
the  inplane  stress  resultants  N  =  {iVn  N22  N12},  bending  moments  M  =  {Mn  M22  Mn}, 
and  transverse  shear  stress  resultant  Q  =  {Q13  Q23},  for  each  region  of  the  multi-plate  model 
as. 
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where 
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GS.  =  /  S„S„Elf(X!,)dx, 

Jti 


In  addition  to  a  beam  element,  BUCKDEL  implements  a  three  noded  triangular  curved 
shell  element.  The  shell  element  is  described  in  the  curvilinear  coordinate  system  x  —  y  and 
the  area  coordinates  are  used  for  field-description.  Accordingly  we  have, 


3 

{x  y  1}  =  Y,U{x  y  \}i  (5) 

2=1 

Inverting  the  above  relationship,  we  get. 


Li  —  {o^nx  +  ai2y  +  Ois)  (6) 

where 

yj  J/fc}  Ui2  —  ^jt  Ot3  =  Xjyk  ^kVj 

^  ^  {^2yz  -  Xzy2  +  Xzyi  -  xiyz  •+■  0:1^2  -  X2yi) 

and  j  =  2, 3,1;  A:  =  3, 1,2  as  i  =  1,2,3. 

The  inplane  displacements  and  the  transverse  shear  strains  need  to  satisfy  C^-continuity 
while  the  transverse  deflection  need  to  satisfy  C^-continuity  in  the  present  formulation.  The 
independent  field  variables  tt,  v,  w,  7x2  and  jyz  are  expressed  in  terms  of  the  nodal  degrees 
of  freedom  Ui,  Vi,  Wi,  =  {-w,y  )i,  Qy.  =  {w,x)i,  Jxzi  and  jy^.  as. 
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where 

Nu  =  Li  +  ^Lj  +  LjLk-LiL^j-LiLl 

N2i  =  dki  {l^Lj  -H  -LiLjL^  +  Gji  (jijLk  +  -LiLjLk^ 

Nzi  =  ak2  {l^Lj  -1-  -LiLjL/^  +  aj2  (j^Lk  +  '^LiLjL^  (8) 

are  the  cubic  polynomials  for  the  transverse  deflection. 

In  the  above  element  formulations,  the  inter-element  C®— continuity  is  exactly  satisfied 
for  all  the  field  variables.  However,  the  inter-element  (7^— continuity  required  for  the  trans¬ 
verse  deflection,  in  case  of  shell  element,  is  satisfied  a  'posteriori  in  a  weak  form  using  the 
Hu-Washizu  variational  principle. 
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BUCKDEL  uses  an  automated  method  to  follow  the  post-buckling  paths  of  the  damaged 
structure.  Automated  post-buckling  solution  involves:  detection  of  possible  instability  in 
solution  and  elimination  of  possible  path-retracing;  classification  of  the  detected  instabilities; 
and  computation  of  the  post-through  buckling  solution(s).  Solution  instabilities  are  detected 
by  monitoring  the  rank  of  the  tangent  stiffness  matrix.  Whenever  the  determinant  of  the 
tangent  stiffness  matrix  changes  its  sign  the  solution  senses  possible  instabilities  in  that  range 
of  load  and  changes  the  sign  of  the  next  load  increment  to  avoid  path-retracing.  Through 
a  cycle  of  iterations,  location  of  instabilities  are  identified  as  the  load  levels  for  which  the 
tangent  stiffness  becomes  singular.  The  tangent  stiffness  is  often  scaled  to  minimize  numerical 
errors.  The  identified  instability  points  are  then  classified  as  limit  points  or  bifurcation  points 
using  some  simple  and  cost-effective  rules  [Huang  and  Atluri  (1995)].  If  the  instability  point 
is  a  limit  point,  the  arc-length  continuation  is  enough  to  obtain  post-buckling  solution  path. 
However,  if  the  instability  point  happens  to  be  a  bifurcation  point,  the  strategies  described 
in  detail  in  Huang  and  Atluri  (1995)  are  used  to  trace  the  appropriate  post-buckling  solution 
branch.  The  nonlinear  fundamental  state  between  the  two  solution  points  n  -  1  and  n  in 
the  neighborhood  of  the  bifurcation  point  is  linearized  to  obtain  the  asymptotic  solution 
for  obtaining  an  approximate  critical  buckling  load  factor.  A  linear  combination  of  the 
normalized  eigen-vector  associated  with  the  critical  buckling  load  factor  and  its  orthogonal 
counterpart  is  used  to  used  to  determine  the  initial  post-buckling  paths. 

3.2.5  Look  Up  Tables 

The  finite  element  based  methodologies  for  modeling  the  various  types  of  damage  being 
considered  in  this  project  are  quite  general  (i.e.  they  are  applicable  to  complex  structural 
geometries  and  metallic  and  composite  materials).  They  are  also  computationally  efficient 
and  very  accurate.  Nevertheless,  it  is  sometimes  useful  to  have  access  to  cataloged  solutions 
(i.e.  stress  intensity  factors,  energy  release  rates,  etc.),  which  can  be  used  to  obtain  estimates 
of  the  residual  strength  and  life  of  common  structural  geometries.  In  addition  to  saving 
computational  time,  these  ‘look  up’  tables  can  be  used  to  provide  sanity  checks  on  the  finite 
element  calculations.  There  are  existing  programs  which  are  essentially  ‘look  up’  tables 
of  stress  intensity  factors  and  material  property  data.  Two  such  programs  are  AFGROW 
(formally  MODGRO)  and  NASGRO.  It  is  not  the  intent  of  the  present  project  to  duplicate 
effort  already  expended  on  programs  such  as  AFGROW  and  NASGRO.  Rather,  some  of  the 
more  common  damage  scenarios  such  as  the  center  cracked  panel,  corner  crack  at  a  hole,  etc. 
will  be  added  to  the  damage  tolerance  module.  By  doing  so,  the  user  of  the  damage  tolerance 
module  will  be  able  to  access  parameters  such  as  stress  intensity  factors  for  these  common 
damage  scenarios  without  having  to  go  through  a  finite  element  analysis.  Besides  stress 
intensity  factor  solutions  for  structures  containing  cracks,  buckling  loads  for  structures  with 
holes  of  various  sizes  and  shapes  will  be  available  to  the  user  through  the  ‘look  up’  tables. 
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3.2.6  Fatigue  Crack  Growth 


The  problem  of  fatigue  crack  growth  is  of  considerable  practical  importance  when  designing 
a  structure  which  satisfies  the  DTR.  To  successfully  employ  damage  tolerance  principles,  an 
accurate  determination  is  required  of  the  number  of  load  cycles  to  failure  in  a  component. 
These  estimates  will  have  a  significant  influence  in  the  design  and  maintenance  of  a  safe 
structure  such  as  in  the  scheduling  of  inspection  intervals. 

Fatigue  crack  growth  frequently  occurs  when  a  flawed  component  is  subject  to  some 
form  of  cyclic  loading.  Here  the  crack  growth  is  termed  subcritical  since,  due  to  the  cyclic 
loading,  it  takes  place  at  stress  intensity  factor  levels  that  are  less  than  the  fracture  toughness 
of  the  material.  The  form  of  the  cyclic  loading  is  also  of  great  importance  such  as  whether  it 
has  constant  amplitude  or  has  a  variable  amplitude.  The  damage  tolerance  module  will  have 
the  capability  to  model  fatigue  crack  under  conditions  of  constant  amplitude  and  variable 
amplitude  loading.  This  capability  will  be  limited  to  self  similar  (i.e.  Mode  I)  crack  growth. 

The  crack  growth  calculations  will  be  performed  based  on  stress  intensity  factors  ob¬ 
tained  either  by  the  FEAM  or  from  the  ‘look  up’  tables  containing  cataloged  solutions.  For 
reasons  of  computational  efficiency,  the  ‘look  up’  table  method  is  preferred  if  a  solution  is 
available.  For  geometries  not  having  cataloged  solutions,  use  will  be  made  of  the  Finite 
Element  Alternating  Method.  The  FEAM  has  made  the  use  of  finite  elements  in  fatigue 
calculations  feasible.  The  reason  being  that  the  global  stiffness  matrix  of  the  finite  element 
model  is  assembled  and  decomposed  only  once,  since  the  stiffness  of  the  uncracked  structure 
remains  the  same  for  all  crack  sizes.  In  other  finite  element  approaches,  such  as  those  us¬ 
ing  singular/hybrid  type  special  crack-tip  finite  elements  or  using  EDI  methods,  the  cracks 
must  be  modeled  explicitly.  Therefore,  the  global  stiffness  matrix  must  be  assembled  and 
decomposed  for  each  crack  size.  The  assembly  and  decomposition  of  the  global  stiflhiess 
matrix  accounts  for  about  80%  of  the  computational  time  required  in  a  typical  finite  element 
analysis.  Given  that  the  FEAM  has  to  perform  this  operation  only  once  during  a  fatigue 
analysis,  the  benefits  of  this  approach  over  other  finite  element  techniques  for  fatigue  crack 
growth  are  readily  apparent. 

Load  spectra,  in  terms  of  ASTROS  load  cases,  will  be  provided  to  the  damage  toler¬ 
ance  module  by  the  USAGE  module  described  in  Nees  (1995).  Since  the  USAGE  module 
defines  aircraft  maneuvers  in  terms  of  ASTROS  load  cases,  the  load  spectra  at  any  point 
in  the  aircraft  structure  can  be  easily  determined  by  extracting  stresses  at  that  point  from 
the  ASTROS  solution  database.  Use  will  be  made  of  USAGE  features  such  as  repeating 
and  blocking  of  data  to  minimize  storage  and/or  computatational  requirements  of  variable 
amplitude  crack  growth  calculations. 

Numerous  studies  have  been  conducted  on  the  characteristics  of  fatigue  crack  growth. 
It  has  been  established  that  when  the  plastic  or  inelastic  zone  in  the  vicinity  of  the  crack  is 
small,  then  the  stress  intensity  factor  is  the  governing  parameter  during  crack  growth  [Paris 
and  Erdogan  (1963)].  In  general,  the  crack  growth  rate  is  a  function  of  stress  intensity  factor 
change,  which  is  given  by: 

=  Kfnax  ^min  (9) 
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where  Kmax  is  the  maximum  stress  intensity  factor  during  the  load  cycle  and  Kmin  is  the 
minimum  value  of  the  stress  intensity  factor.  Based  on  the  minimum  and  maximum  loads, 
it  is  customary  to  define  the  parameter  R  where: 


—  ^min/ ^max  (Ifi) 

An  important  point  must  be  made  in  relation  to  the  functional  relationship  between  the 
fatigue  crack  growth  rate  and  the  change  in  stress  intensity  factor.  This  crack  growth  rate 
function  can  be  partitioned  into  three  separate  regions.  At  low  values  of  AK,  there  is  very 
little  crack  growth  with  a  negligible  crack  growth  rate.  Therefore,  it  can  be  stated  that 
there  exists  a  AK  below  which  there  is  no  crack  growth.  This  quantity  is  referred  to  as 
the  threshold  stress  intensity  factor  and  is  denoted  as  {AK)th.  At  higher  values  of  AK, 
crack  growth  takes  place.  It  has  been  observed  experimentally  that  this  curve,  relating  the 
crack  growth  rate,  da/dN,  to  AK,  is  usually  linear  on  a  log-log  plot  and  this  corresponds 
to  a  power  law  relation  between  da/dN  and  AK.  This  is  commonly  referred  to  as  the  Paris 
relation  and  is  given  as  [Paris,  Gomez  and  Anderson  (1961)]: 

^=C(AKr  (11) 

where  a  is  the  crack  length  and  N  is  the  number  of  load  cycles.  The  quantities  C  and 
n  are  material  dependent  constants.  At  higher  values  of  AK  the  stress  intensity  factor  is 
approaching  the  fracture  toughness  of  the  material,  Kc-  The  crack  growth  rate  will  increase 
significantly,  eventually  leading  to  the  onset  of  rapid  unstable  crack  growth 

Recognizing  that  several  distinct  phases  in  fatigue  crack  growth  exist,  a  more  general 
form  for  the  relationship  between  crack  growth  rate  and  stress  intensity  factor  is  expressed 
as: 

da  C{1  -  RrjAKnAK  - 

dN  {(1  -  R)Kc  -  AK}<3  ^  ’ 

where  m,p  and  q  are  constants  that  relate  to  the  particular  crack  growth  relation  that  is 
being  used.  By  assigning  different  values  to  these  quantities  some  of  the  well  known  crack 
growth  relations  can  be  recovered.  For  example,  when  m  =  p  =  q  =  0,  the  Paris  relation 
is  obtained.  The  Forman  relation  [Forman,  Kearney  and  Engle  (1967)],  which  accounts  for 
high  crack  growth  rates  and  instability,  is  recovered  by  setting  m  =  p  =  0  and  q  =  1.  By 
setting  p  =  g  =  0  and  m  =  (M^  -  l)n,  the  Walker  relation  [Walker  (1970)]  is  derived,  where 
Miy  is  an  exponent  in  the  Walker  relation. 

With  stress  intensity  factor  solutions  and  a  crack  growth  relation  in  hand,  the  ultimate 
objective  of  fatigue  crack  growth  calculations,  to  calculate  the  number  of  cycles  for  a  crack 
to  grow  by  a  specified  amount,  can  be  carried  out.  This  is  done  by  integrating  the  growth 
relation  [Equation  (11)  or  (12)].  For  example,  equation  (11)  is  integrated  to  give: 


1  /•“!  da 

AJao  Ai^ 


(13) 


where  Cq  and  Oi  are  the  initial  and  final  crack  lengths.  This  integration  is  done  numerically  by 
sub-dividing  the  crack  growth  distance  into  a  number  of  intervals  and  calculating  the  stress 
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intensity  factor  at  each  interval.  A  trapezoidal  rule  is  then  used  to  calculate  the  number  of 
cycles.  All  of  these  steps  will  be  carried  out  automatically  by  the  damage  tolerance  module. 

Self-similar  growth  of  through  skin  cracks  is  straightforward  to  carry  out  using  the 
above  procedure.  However,  for  an  elliptical  or  part  elliptical  crack  since  the  stress  intensity 
factor  distribution  is  generally  not  constant  along  the  elliptical  crack  front,  the  crack  growth 
rates  will  not  be  the  same  in  every  direction.  Consequently,  the  shape  of  the  crack  will 
change.  It  has  been  often  observed  in  practice  that  a  semi-elliptical  surface  crack,  initially 
having  a  small  aspect  ratio,  will  have  a  larger  crack  growth  rate  in  the  minor  axis  direction. 
Therefore,  in  the  damage  tolerace  module  crack  growth  will  be  based  on  the  extension  of 
the  minor  axis  by  a  specified  amount.  The  number  of  cycles  for  this  is  then  calculated  be 
evaluating  the  stress  intensity  factor  at  each  sub-interval  as  described  earlier.  The  procedure 
to  adjust  the  major  axis  length  is  as  follows.  Based  on  the  extension  of  the  minor  axis  in  a 
given  interval,  an  estimate  is  made  of  the  extension  of  the  major  axis.  The  cycles  for  crack 
growth  in  these  two  directions  are  calculated  using  Equation  (13).  Since  these  are  different, 
in  general,  the  crack  extension  in  the  major  axis  direction  is  adjusted  such  that  the  number 
of  cycles  is  the  same  for  both  directions  in  a  given  interval.  This  is  done  by  simple  linear 
interpolation /extrapolation. 

3.3  Integration  of  the  Damage  Tolerance  Module  with  ASTROS 

This  section  will  address  the  integration  of  the  damage  tolerance  module  into  ASTROS.  An 
important  point  which  must  be  addressed  is  the  level  of  integration  between  the  damage 
tolerance  module  and  ASTROS.  Here,  level  of  integration  refers  to  the  extent  to  which  the 
damage  tolerance  module  will  make  use  of  existing  ASTROS  modules  which  handle  tasks  such 
as  I/O,  memory  management,  large  matrix  operations,  etc.  This  subject  was  investigated 
during  Phase  I  and  it  was  concluded  that  the  damage  tolerance  module  will  be  designed 
to  run  in  a  ‘stand  alone’  mode  or  as  an  analysis  module  in  ASTROS.  This  conclusion  was 
reached  for  the  following  three  reasons. 

1.  From  a  development  standpoint,  complete  control  over  the  source  code  is  needed.  If 
the  software  can  be  compiled  and  executed  independent  of  ASTROS,  this  will  assist  in 
debugging,  verification,  and  maximization  of  computational  efficiency.  Further,  future 
porting  and  maintenance  of  the  software  will  be  much  simpler. 

2.  For  the  purposes  of  commercialization  of  the  software  (which  is  a  requirement  of  Phase 
II),  some  potential  customers  may  need  only  the  damage  tolerance  module  and  not  the 
global  modeling  and  multidisciplinary  design  and  optimization  capabilities  of  ASTROS. 
Thus,  they  may  not  be  willing  to  pay  the  extra  costs  associated  with  ASTROS  just  to 
be  able  to  make  use  of  the  damage  tolerance  module. 

3.  Nothing  is  lost  by  developing  the  damage  tolerance  module  to  run  in  both  a  ‘stand 
alone’  mode  and  as  an  ASTROS  analysis  module  as  opposed  to  being  developed  to  run 
only  with  ASTROS. 
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The  software  developed  in  this  SBIR  project  will  consist  of  four  components.  These  are 
ASTROS,  the  USAGE  module,  an  interface  module  and  the  damage  tolerance  module.  The 
software  architecture  is  shown  schematically  in  Fig.  15.  The  ASTROS  system  will  be  rebuilt 


Figure  15:  Software  Architecture 

with  a  modified  MAPOL  sequence,  new  relational  entities,  new  bulk  data  entities,  new  error 
messages  and  the  new  module  definition  so  that  it  will  call  the  damage  tolerance  module 
when  instructed  to  do  so  by  the  global  model  input.  The  USAGE  module,  described  in  Nees 
(1995),  will  be  used  to  generate  load  spectra  in  terms  of  ASTROS  load  cases.  The  interface 
and  damage  tolerance  module  will  be  coded  in  FORTRAN  77  by  Knowledge  Systems.  The 
interface  module  will  facilitate  the  exchange  of  data  between  ASTROS,  the  USAGE  module 
and  the  damage  tolerance  module.  The  data  flowing  from  ASTROS  to  the  damage  tolerance 
module  will  be  information  on  the  contstraints  and  associated  sensitivities  to  be  evaluated 
for  each  local  model,  current  values  of  the  design  variables,  and  boundary  conditions  to  be 
applied  to  the  local  models.  The  data  flowing  from  the  damage  tolerance  module  to  ASTROS 
will  be  the  constraints  and  sensitivities  which  were  evaluated  by  the  damage  module. 
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Panel  A 


Figure  16:  Global  Finite  Element  Model  of  the  Intermediate  Complexity  Wing  (ICW) 

3.3.1  Illustration:  BUCKDEL  Integration  with  ASTROS 

As  an  explicit  illustration  of  how  this  integrated  system  will  work,  this  section  describes 
how  the  BUCKDEL  software  [Pipkins  and  Atluri  (1996)],  developed  in  Phase  I,  would  be 
used  to  introduce  panel  buckling  constraints,  in  the  presence  of  delaminations,  to  the  design 
of  the  Intermediate  Complexity  Wing  (ICW).  In  this  illustration,  BUCKDEL  assumes  the 
role  of  the  damage  tolerance  module  of  Fig.  15.  The  ICW  model  is  described  in  Sections 
4.7  and  4.8  of  the  ASTROS  Applications  Manual.  Fig.  16  shows  the  global  finite  element 
model  of  the  ICW.  The  model  consists  of  39  rod  elements,  55  shear  panel  elements,  62 
quadrilateral  plate  elements  and  2  triangular  plate  elements.  The  substructure  material  is 
modeled  as  aluminum  while  the  wing  skins  are  made  of  a  graphite/epoxy  composite.  The 
ICW  examples  in  the  ASTROS  Applications  Manual  are  designing  a  minimum  weight  wing 
under  strength  and  flutter  constraints.  In  this  illustration,  panel  buckling  constraints  in  the 
presence  of  delaminations  are  also  introduced  into  the  design.  The  panel  buckling  constraint 
is  applied  to  Panel  A  (see  Fig.  16).  A  local  finite  element  model  of  Panel  A  (Fig.  17)  is 
developed  and  will  be  read  as  input  by  BUCKDEL,  when  it  is  called  by  ASTROS. 

The  integration  of  BUCKDEL  into  version  11  of  ASTROS  is  quite  easy  due  to  the 
presence  of  an  existing  panel  buckling  constraint  for  rectangular  panels  with  no  damage. 
With  just  a  slight  modification  to  the  MAPOL  sequence,  BUCKDEL  can  be  used  in  place  of 
this  panel  buckling  constraint  capability  which  is  included  in  version  11  of  ASTROS.  First, 
the  Module  Definitions  file  (MODDEF.DAT)  is  modified.  By  doing  so,  BUCKDEL  can  then 
be  called  by  the  MAPOL  code.  Next  the  MAPOL  sequence  is  edited.  MAPOL  calls  to 
the  engineering  application  modules  PBKLEVAL  and  PBKLSENS  are  replaced  with  calls 
to  BUCKDEL.  The  arguments  passed  to  and  from  BUCKDEL  would  be  the  same  as  those 
passed  to  and  from  PBKLEVAL  and  PBKLSENS  except  that  a  flag  to  inform  BUCKDEL 
whether  it  is  to  evaluate  a  constraint  or  constraint  sensitivity  is  also  passed.  With  the 
above  changes  made  to  the  MODDEF.DAT  file  and  the  MAPOL  sequence,  the  SYSGEN 
program  is  run.  SYSGEN  will  generate  a  new  XQDRIV  subroutine,  which  provides  the 
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Figure  17:  Local  Finite  Element  Model  of  Panel  A  with  Delamination 

actual  link  between  the  MAPOL  call  and  the  FORTRAN  routines  that  make  up  BUCKDEL. 
The  new  XQDRIV  subroutine  must  be  compiled  and  linked  with  the  ASTROS  object  file  to 
produce  the  new  executable  version  of  ASTROS  which  can  make  use  of  the  damage  modeling 
capabilities  in  BUCKDEL. 

From  this  point,  the  use  of  BUCKDEL  to  evaluate  the  panel  buckling  constraint  and 
associated  sensitivities  for  panel  A  of  the  ICW  model,  with  account  taken  for  the  presence 
of  a  delamination,  is  exactly  the  same  as  the  use  of  the  panel  buckling  constraint  capability 
shipped  with  version  11  of  ASTROS.  Namely,  the  user  must  specify  a  DCONBK  bulk  data 
card  for  panel  A.  After  reading  the  DCONBK  bulk  data  card,  ASTROS  will  call  BUCKDEL 
at  the  appropriate  points  in  the  MAPOL  sequence.  When  called,  BUCKDEL  will  read  the 
local  finite  element  model  data  for  panel  A,  extract  necessary  information  such  as  panel  loads 
from  the  ASTROS  database  and  carry  out  the  constraint  or  constraint  sensitivity  evaluation. 
After  evaluating  this  information,  BUCKDEL  will  store  this  information  in  the  ASTROS 
database  for  use  by  the  optimization  routine. 
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Appendix  A:  Schwartz-Neumann  Alternating  Method 

The  Schwartz-Neumann  alternating  method  is  based  on  the  superposition  principle.  The 
solution  on  a  given  domain  is  the  sum  of  the  solutions  on  two  other  overlapping  domains, 
with  part  of  the  boundary  conditions  as  unknowns.  The  alternating  method  can  be  viewed 
as  the  fixed  point  iteration  scheme  used  to  solve  these  unknown  boundary  conditions.  Based 
on  this  point  of  view,  we  can  perform  a  convergence  study.  The  alternating  method  con¬ 
verges  unconditionally  when  there  are  only  traction  boundary  conditions  specified  on  the 
body.  The  convergence  criterion  for  mixed  boundary  value  problems,  where  there  are  ap¬ 
plied  displacement  boundary  conditions  as  well  as  traction  boundary  conditions,  is  discussed 
in  the  following.  Compare  the  work  done  by  the  applied  forces  in  the  following  two  cases. 
In  the  first  case,  arbitrary  displacement  conditions  exist  on  the  surfaces  of  the  cracks  in  the 
cracked  finite  body,  while  all  the  boundary  conditions  elsewhere  are  replaced  by  homogeneous 
boundary  conditions,  i.e.  remove  all  tractions  and  reduce  all  the  applied  displacements  to 
zero  magnitude.  In  the  second  case,  the  same  displacement  conditions  exist  on  the  surfaces 
of  the  cracks  in  the  infinite  domain.  If  the  work  done  in  the  cracked  finite  body  is  always 
smaller  than  twice  of  the  work  done  in  the  infinite  domain,  the  alternating  method  converges. 
Otherwise,  it  does  not.  For  most  practical  problems,  this  ratio  is  close  to  one.  Thus,  the 
alternating  method  converges  rapidly,  as  discussed  in  detail  in  the  following  section. 

A.l  Superposition  Principle  and  the  Alternating  Method 

Consider  n  cracks  in  a  body  of  a  finite  size.  The  crack  surfaces  which  are  traction  free,  are 
denoted  collectively  as  Fc.  Let  the  boundary  of  the  finite  domain(not  including  the  crack 
surface)  be  T,  of  which  the  boundary  with  prescribed  tractions  t°  is  T*,  and  the  boundary 
with  prescribed  displacements  w®  is  r„.  It  is  clear  that  T  =  r„  U  Ft. 

The  alternating  method  uses  the  following  two  simpler  problems  to  solve  the  original 
one.  The  first  one,  denoted  as  Pan  a  [shown^  in  Fig.  18(c)],  is  that  of  the  same  n  cracks  in  the 
infinite  domain  subjected  to  the  unknown  crack  surface  loading  T.  The  second  one,  denoted 
as  PpEM  [shown  in  Fig.  18(b)],  has  the  same  finite  geometry  as  in  the  original  problem  except 
that  the  cracks  are  ignored.  The  boundary  F^  of  Pfem  has  the  prescribed  displacement  w, 
while  the  boundary  Ft  has  the  prescribed  traction  t.  The  prescribed  displacements  and 
tractions  are  different  from  those  in  the  original  problem  in  general.  Because  of  the  absence 
of  the  cracks,  the  problem  Pfem  can  be  solved  much  easier  by  the  finite  element  method  (or 
the  boundary  element  method). 

To  solve  the  original  problem,  Pqrg  (shown  in  Fig.  18(a)),  the  crack  surface  loading  T, 
the  prescribed  displacement  u  and  the  traction  t  must  be  found  such  that  the  superposition 
of  the  two  alternative  problems  Pan  a  and  Pfem  yields  the  original  one,  Pqrg-  The  detailed 
procedures  to  find  these  boundary  conditions  are  described  as  the  following. 

In  the  uncracked  body  problem  Pfem,  the  tractions  T  at  the  location  of  the  cracks  in 
^Fig.  18  only  illustrates  one  crack.  Many  cracks  may  be  present. 
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Figure  18:  Superposition  Principle  for  Finite  Element  Alternating  Method 

the  cracked  body  Porg  can  be  solved,  for  any  given  boundary  loads  u  and  t,  using  the  finite 
element  method.  Due  to  the  linearity  of  the  problem,  the  solution  can  be  denoted  as 

T  =  +  KH  (14) 

where  and  are  linear  operators. 

Similarly,  the  tractions  on  boundary  F*  and  the  displacements  on  boundary  r„ 
can  be  found  in  the  infinite  domain  Pan  a  for  the  given  crack  surface  load  T,  which  is  the 
same  as  the  crack  surface  traction  obtained  in  the  Pfem-  The  solution  can  be  denoted  as 

=  H^T  (15) 

(16) 

where  and  are  also  linear  operators.  Subtract  the  solution  for  Pana  from  the  one 

for  Pfem-  The  resulting  solution  has  zero  tractions  at  the  location  of  the  crack  surfaces.  To 
ensure  that  the  resulting  solution  has  the  same  boundary  conditions  on  F,  the  Eq.  (17)  and 
Eq.  (18)  must  be  satisfied. 

u  =  n"  +  (17) 

t  =  t°  +  e  (18) 

The  unknown  tractions  T,  t  and  unknown  displacement  u  can  be  solved  using  these 
equations  [Eq.  (14)  through  Eq.  (18)].  Eliminate  w,  and  t,  by  substituting  Eq.  (17), 
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Eq.  (18),  Eq.  (15)  and  Eq.  (16)  into  Eq.  (14)  to  obtain  the  following  equation  for  the  traction 
T.  _ 

I  -  T  =  +  KH°)  (19) 

Eliminate  and  T  to  obtain  the  following  equation  for  the  unknown  traction  t  and 

unknown  displacement  u. 

“}  =  {“•}  (20) 

where  _  _ 

_  r  ■ 

~  [  K^K'^  K*K^ 

and  I  is  the  identity  operator. 

Similarly,  we  can  obtain  the  following  linear  system  for  traction  and  displacement 

(I-A)X  =  y  (21) 


where 


Y  =  A 


It  is  possible  to  solve  these  equations  directly  to  obtain  the  tractions  T,  t  and  displace¬ 
ment  u.  But  this  involves  the  evaluation  of  and  K*,  which  requires  solving  the  traction 
T  at  the  location  of  the  uncracked  body  subjected  to  all  different  loading  patterns  u  and  t. 
We  have  to  solve  the  uncracked  body  problem  a  larger  number  of  times,  of  the  same  order 
as  that  of  the  total  number  of  degrees  of  freedom  of  the  boundary  nodes,  using  the  finite 
element  method.  Thus,  it  can  be  very  expensive  to  find  X  by  solving  directly  the  linear 
system  (I  —  A)X  =  Y.  A  fixed  point  iteration  scheme  can  be  used  to  solve  this  linear 
system.  The  iterative  scheme  can  be  devised  as: 


z  =  0,1,2,...,  oo 


where  =  If  this  procedure  converges,  the  solution  is 


(22) 


i=l 


Since  _ 

){jr“.  K‘} 

the  iterative  scheme  Eq.  (22)  is  equivalent  to  the  following  alternating  scheme 

-t-  (23) 
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Figure  19:  Subtract  \Pfem  from  Pana  to  Obtain  the  Solution  for  Pres  which  has  Homo¬ 
geneous  Boundary  Condition  on  F 


for  i  =  0, 1, 2, ,  00.  In  this  case,  the  uncracked  body  problem  is  solved  only  a  few  times, 
because  this  fixed  point  iteration  scheme  converges  quickly  for  practical  problems.  Therefore, 
the  alternating  method  is  much  more  efficient  than  solving  the  linear  system  directly.  But  it 
should  be  noticed  that  it  may  not  be  necessary  to  use  the  alternating  method  in  some  cases. 
It  can  be  more  efficient  and  accurate  to  solve  directly  when  multiple  crack  solutions  are 
constructed  from  that  for  a  single  crack.  This  will  be  discussed  in  detail  in  a  later  section. 


A.2  Convergence  of  the  Alternating  Method 

First  it  is  shown  that  I  -  A  is  not  singular  and  the  linear  system  Eq.  (21)  has  a  unique 
solution.  Suppose  J  —  A  is  singular.  Then,  there  must  exist  a  non-zero  X  such  that 
{I  —  A)X  =  0,  which  means  that  there  exist  non-zero  and  and  therefore  a  non-zero 
T,  such  that 

T  =  -i-  KH^ 
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In  this  case  the  analytical  solution  and  the  finite  element  solution  have  the  same  displacement 
u“  on  boundary  r«  and  the  same  traction  on  boundary  Tf.  Subtracting  the  analytical 
solution  from  the  FEM  solution,  we  obtain  the  solution  for  the  following  problem.  The  entire 
boundary  F  is  free  of  external  loadings  as  well  as  the  crack  surfaces.  But,  the  FEM  solution 
gives  zero  displacements  for  the  crack  surfaces,  while  the  analytical  solution  gives  non-zero 
displacements  for  the  crack  surfaces  because  of  the  non-zero  T.  Thus,  the  resulting  solution 
has  non-zero  displacements  at  the  crack  surfaces.  This  is  a  contradiction  because  the  cracks 
cannot  be  opened  without  any  external  load.  Consequently,  I  —  Ais  not  singular. 

The  fixed  point  iteration  scheme  Eq.  (22)  converges  if  all  the  eigenvalues  of  A  are  in 
the  open  interval  (—1,1).  The  scheme  of  Eq.  (22)  converges  since  the  eigenvalues  of  A  are 
in  (-1, 1)  for  most  problems  of  practical  interest.  The  eigenvalues  of  A  are  smaller  than  1. 
Let  X A  be  an  eigenvector  of  A  corresponding  to  the  eigenvalue  A. 

T  =  K'^(u^)  +  K\t^) 

Xu^  =  WT 
Xt^  =  K*T 

The  solution  Prbs^  shown  in  Fig.  19(a),  is  obtained  by  subtracting  A  times  the  FEM  solu- 
tion(Fig.  19(c))  from  the  analytical  solution(Fig.  19(b)).  Here,  u  =  0  and  t  =  0  on  F  and 
the  crack  surface  loading  is  (1  —  A)r,  while  the  displacements  at  the  crack  surfaces  are  the 
same  as  those  in  the  analytical  solution.  If  the  work  done  in  opening  the  cracks  in  the  infinite 
domain  is  W,  the  work  done  in  opening  the  cracks  in  the  finite  domain  (with  the  boundary 
condition  u  =  0  and  t  =  0)  is  (1  —  X)W,  which  is  equal  to  the  strain  energy  stored  in  the 
body.  It  must  be  positive.  Thus,  A  <  1. 

It  can  be  shown  that  A  >  0  in  the  absence  of  the  prescribed  displacement  boundary 
conditions.  In  this  case,  the  resulting  solution  from  the  subtraction  has  zero  tractions  at  the 
boundary  F .  Apply  additional  load  At  to  the  boundary  F  with  the  crack  surfaces  fixed.  The 
stress  state  in  the  body  will  be  the  same  as  that  in  the  analytical  solution  described  above, 
after  this  additional  loading  is  applied.  This  procedure  of  adding  load  on  the  boundary  F  is 
exactly  the  same  as  that  in  the  FEM  solution  for  the  uncracked  body,  except  that  the  load 
level  is  A  times  of  that  in  the  FEM  solution,  because  the  crack  surfaces  are  fixed.  Therefore, 
the  work  done  by  the  additional  load  is  positive.  Consequently,  (1  -  X)W  <  W  and  A  >  0. 
So,  the  alternating  method  converges  for  cracks  in  finite  domains  with  arbitrary  shapes  and 
arbitrary  traction  boundary  conditions. 

In  general,  the  eigenvalue  A  can  be  smaller  than  zero  for  mixed  boundary  problems.  It 
is  greater  than  -1  only  if  (1  —  A)IF  <  2W.  Thus,  the  convergence  criterion  for  the  alternating 
method  for  the  general  case  with  mixed  boundary  conditions  can  be  stated  as  follows.  The 
alternating  method  [Eq.  (22)]  converges  if  the  crack  surface  loads  do  less  work  in  the  finite 
domain,  with  the  homogeneous  boundary  condition  «  =  0  and  t  =  0  on  F,  than  twice 
as  much  as  they  do  in  the  infinite  domain  for  any  arbitrary  distribution  of  crack  surface 
displacements  (see  Fig.  20). 

Quick  convergence  can  be  expected  for  most  of  the  practical  applications.  For  any  crack 
surface  displacements,  the  displacements  and  stresses  at  a  point  decay  rapidly  as  the  point 
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FEAM  converges 
if  the  work  done  in  PpjisjrrE 
less  than  twice  of  the  work  done  in 


the  infinite  body  with  any 
arbitrary  crack  surface 
displacement 


the  finite  body  with  the 
same  crack  surface 
displacement 


(a) 


(b)  PpiNITE 


Figure  20:  Convergence  Criterion 

moves  away  from  the  cracks.  Thus,  the  work  done  in  the  finite  domain  with  the  homogeneous 
boundary  condition  is  very  close  to  the  work  done  in  the  infinite  domain.  This  implies  that 
the  eigenvalues  of  A  are  very  small  and  the  fixed  point  iteration  converges  rapidly.  Indeed, 
all  mixed  boundary  value  problems  we  have  solved  (for  both  2D  and  3D  problems)  to  date 
using  finite  element  alternating  method  have  converged. 


A. 3  Summary  of  FEAM  Procedure 

The  alternating  procedure  defined  in  Eq.  (22)  can  be  translated  into  the  following  simple 
procedure.  Refer  to  Fig.  18. 

1.  Solve  PpEM  with  the  given  load  on  the  boundary  F.  Solve  for  the  tractions,  which  are 
used  to  close  the  cracks.  Denote  the  solution  as  where  1  indicates  that  this  is 

the  solution  for  the  first  iteration. 

gFEM.  =  K'^u°  +  KH° 
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2.  Reverse  the  crack  surface  traction  obtained  in  the  previous  step  and  apply  it  as  the 
load  on  the  crack  surfaces  and  solve  the  Pan  a-  Denote  the  solution  as 


3.  Find  the  tractions  on  the  boundary  F*  and  the  displacements  on  the  boundary  F^  from 
the  analytical  solutions  obtained  in  the  previous  step.  Reverse  them  as  the  load  for 
PpEM-  Find  the  crack  closing  tractions  from  the  solution  . 

gFEM.  2’(2)  = 


4.  Repeat  the  step  2  and  3  until  the  residual  load  is  small  enough  to  be  ignored. 


cFEM, 

• 


T^i+l)  _  4. 


for  z  =  2, 3, . . .. 


The  solution  to  the  original  problem  is  the  summation  of  all  those  obtained  in  the 
alternating  procedure,  i.e. 

S  =  f2{sr^  +  Sf^^)  (25) 

2=1 


A.4  Solutions  for  Multiple  Skin  Cracks 

Solutions  for  multiple  skin  cracks  in  an  infinite  body,  subjected  to  arbitrary  crack  surfaces 
tractions,  can  be  constructed  using  the  solution  for  a  single  crack  in  a  infinite  body  subjected 
to  arbitrary  crack  surface  loading.  Analytical  solutions  for  multiple  embedded  cracks  in  an 
infinite  body  are  available  only  for  some  special  configurations,  such  as  multiple  collinear 
cracks  subjected  to  arbitrary  crack  surface  tractions  [Muskhelishvili  (1953)].  There  are  sev¬ 
eral  implementations  of  finite  element  alternating  method  based  on  such  analytical  solutions 
[Park  and  Atluri  (1993),  Wang  and  Atluri  (1995),  etc.].  It  is  in  general  easier  to  construct 
the  solution  of  multiple  embedded  cracks  in  an  infinite  body  using  the  solution  for  a  single 
embedded  crack.  Solutions  for  arbitrary  distributions  of  cracks  can  be  obtained  using  this 
approach.  It  can  be  more  accurate  and  efficient  to  build  the  multiple  crack  solutions  from 
that  for  a  single  crack  even  when  the  analytical  solution  is  available,  such  as  for  the  multiple 
collinear  cracks  in  an  infinite  domain. 

In  the  context  of  the  finite  element  alternating  method,  it  seems  natural  to  use  the 
Schwartz-Neumann  alternating  method  to  obtain  the  analytical  solution  iteratively.  This 
approach  has  been  used  by  many  authors,  such  as  O’Donoghue,  Nishioka  and  Atluri  (1984) 
and  (1985).  Using  the  alternating  method  and  the  solution  for  the  single  crack  in  the  infinite 
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domain,  residuals  induced  by  closing  the  other  cracks  are  erased  by  reversing  them  and 
applying  them  as  loads  on  the  crack  surfaces.  Consider  n  arbitrarily  distributed  cracks.  Let 

the  given  crack  surface  load  on  k’th  crack  be  T^,  (A:  =  1, 2, _ n).  The  alternating  procedure 

is  outlined  as  the  following. 

1.  Consider  a  single  crack,  located  at  the  same  position  as  that  of  the  fe’th  crack  in  the 
original  problem,  in  the  infinite  domain.  Apply  the  load  on  the  crack  surface. 
Denote  the  analytical  solution  for  this  loading  as 

2.  Compute  the  cohesive  traction  that  is  used  to  close  the  fth.  crack  in  the  original 
problem,  using  the  solution 

■  ^fk  =  j  =  1,2,..., n  and  j 

where  =  T^,  and  [j  =  1,2,  ...,n  and  j  ^  k)  are  linear  operators.  The 
superscript  indicates  that  the  crack  in  the  single-crack  solution  is  at  the  same  location 
as  the  A;’th  crack  in  the  original  multiple-crack  problem. 

3.  Sum  the  crack  closure  tractions  due  to  all  solution  Sj^\  {k  =  1,2, ...  ,n)  to  find  the 
residual  traction.  Reverse  the  residual  tractions  as  loads. 

r<'>  =  -  t  fS>  i  =  l,2,...,n 

k — Ijk^j 

4.  Apply  the  load  on  the  crack  surface  as  in  the  first  step.  Denote  the  solution  as 

Repeat  the  process  of  finding  residual  tractions  until  the  residual  tractions  are 
small  enough  to  be  ignored. 

Tp'^  =  -  t  ;  =  l,2....,n 

k — l,k^j 

for  i  =  1,2, ... . 

The  solution  to  the  multiple  cracks  subjected  to  the  given  load  Tl  (k  =  1,2, . . .  ,n)  is 
the  sum  of  all  the  solutions  involved  in  the  alternating  procedure,  i.e. 

oo  n 

i=l  A;=l 

However,  the  solution  can  be  obtained  using  a  non-iterative  approach  in  a  simpler  and 
more  efficient  fashion.  As  we  have  shown  in  the  previous  section,  an  alternating  method  is 
essentially  a  fixed  point  iteration  scheme  used  to  solve  a  linear  system.  In  the  analysis  of  a 
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given 


Figure  21:  Superpose  the  Single  Crack  Solutions 


cracked  body,  solving  the  linear  system  directly  needs  crack  closure  tractions  for  all  possible 
boundary  loadings.  Thus  the  uncracked  body  problem  has  to  be  solved  a  large  number  of 
times,  of  the  same  order  as  that  of  the  total  number  of  degrees  of  freedom  of  the  nodes  at 
the  boundary,  for  different  loadings.  Only  a  few  iterations  are  needed  when  finite  element 
alternating  method  is  used.  Thus  the  uncracked  body  is  solved  only  a  few  times. 

On  the  other  hand,  these  tractions  can  be  easily  evaluated  in  the  analysis  of  multiple 
cracks  using  the  analytical  solutions.  Since  the  number  of  degrees  of  freedom  involved  is 
small,  the  linear  system  can  be  solved  directly. 

The  linear  system  for  solving  the  multiple  crack  problem  is  derived  from  the  superposi¬ 
tion  principle.  Consider  the  superposition  of  n  solutions  of  single  cracks  in  the  infinite  body. 
Each  of  these  n  solutions  involves  only  one  crack.  Denote  the  k'th  solution  as  Sk,  where  the 
crack  is  at  the  same  location  as  that  of  the  fc’th  crack  of  the  original  multiple-crack  problem. 
The  crack  surface  traction  T*  for  the  problem  Sk  (A:  =  1, 2, . . . ,  n)  is  to  be  determined  (see 
Fig.  21). 

The  traction  at  the  location  of  the  j’th  crack  in  the  the  problem  Sk  can  be  found  for 
any  load  T*,  i.e. 

^jk  —  jjk  —  1, 2, . . . ,  n.  (26) 

It  is  noticed  that  =  J  (/i:  =  1, 2, . . . ,  n)  are  identity  operators,  because  the  tractions  at 
the  crack  surfaces  are  the  same  as  the  applied  loads. 
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The  superposition  of  the  n  solutions  should  give  back  the  original  problem,  i.e.  the 
tractions  at  the  locations  of  the  crack  surfaces  should  be  the  same  as  the  given  crack  surface 
loads.  Thus,  the  linear  system  to  be  solved  is 

Y,tii  =  '£KfT„  =  T'’  i  =  l,2,...,n.  (27) 

k=l  k-1 

We  can  denote  collectively  the  undetermined  crack  surface  loads,  (A:  =  1, 2, . . . ,  n),  as  T. 
Similarly,  denote  collectively  the  given  loads,  (A:  =  1, 2, . . . ,  n),  as  T°,  Thus,  Eq.  (27)  can 
be  rewritten  as 

KT  =  T°  (28) 

where  iT  is  a  linear  operator.  Once  the  linear  operator  K  is  evaluated  numerically,  we  can 
solve  the  linear  system  directly  instead  of  using  alternating  method. 

Crack  surface  tractions  have  to  be  discretized  by  a  set  of  linearly  independent  basis 
functions,  such  as  polynomial  functions,  Chebyshev  polynomials,  or  certain  piecewise  con¬ 
tinuous  functions,  since  arbitrary  crack  surface  tractions  can  not  be  handled  directly  by 
numerical  methods. 

Let  the  undetermined  load  T  be  approximated  by  N  basis  functions  Bj  (7  =  1, 2, . . . ,  iV). 

(29) 

i=i 

Similarly,  the  given  load  T°  which  can  be  approximated  as  in  the  following. 

N 

(30) 

i—\ 

We  apply  load  on  the  cracks.  Close  all  other  cracks  except  the  single  crack  on  which 
Bi  has  non-zero  values.  Find  the  tractions  at  the  locations  of  the  n  cracks  of  the  original 
problem  (see^  Fig.  22  ),  using  the  analytical  solution  for  a  single  crack  in  an  infinite  domain. 

tj  =  KBj  j  =  l,2,...,N  (31) 

where  iiC  is  a  linear  operator.  The  tractions  tj  {j  =  1, 2, . . . ,  iV)  are  also  approximated  by 
these  basis  functions. 

N 

tj  ~  ^  tjiBi  jf  =  1, 2, . . . ,  AT  (32) 

i=l 

Once  the  magnitude  Uj  are  evaluated,  the  traction  at  the  location  of  A;’th  crack  can  be 
determined  for  a  crack  surface  load  T^.  Since 

t  =  KT 

^Point  loads  (Delta  functions)  are  used  to  illustrate  the  base  functions  in  the  figure.  Only  some  of  the 
loading  cases  are  illustrated  in  the  figure. 
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cohesive  tractions 


Figure  22:  Evaluate  the  Traction  at  the  Locations  of  Cracks  for  Each  Load  in  Terms  of  Unit 
Basis  Functions 

=  f:TiKBj 

J=1 
N  N 

j=l  2=1 

the  linear  system  Eq.  (28)  can  be  approximated  as 

N  N  N 

=  (33) 

j  =  l  2=1  2=1 

which  leads  to  the  following  linear  system  of  equations  for  the  magnitudes  (y  =  1, 2, . . . ,  n). 

N 

Y,tnTi=Tt  i=l,2,...,Ar  (34) 

j=l 

The  alternating  method  essentially  solves  the  same  approximated  linear  system  with  a  fixed 
point  iteration  scheme. 

The  non-iterative  procedure  to  solve  the  multiple  crack  problem  is  summarized  as  the 
following. 
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1.  Apply  loads  in  terms  of  unit  basis  functions  on  one  of  the  cracks,  ignoring  the  other 
cracks.  Use  the  analytical  solution  for  a  single  crack  to  solve  the  tractions  at  the 
locations  of  all  the  other  cracks  (see  Fig.  22). 

2.  Approximate  these  tractions  in  terms  of  the  linear  combination  of  the  basis  functions. 
Find  the  magnitude  of  each  component.  Thus,  the  coefficients  of  the  linear  system  is 
determined  as  in  Eq.  (32). 

3.  Approximate  the  given  crack  surface  load  T*  in  terms  of  the  linear  combination  of  the 
basis  functions.  Find  the  magnitude  of  each  components. 

4.  Solve  the  linear  system  Eq.  (34)  to  obtain  the  the  loads  applied  on  each  single  crack 
solution. 

5.  Superpose  the  n  single  crack  solutions  to  form  a  solution  for  the  original  problem. 

The  coefficients  of  the  linear  system  remain  the  same  in  the  analysis  of  the  same  cracks 
under  different  loadings,  because  they  depend  only  on  the  crack  configuration  and  the  basis 
functions.  Thus,  the  linear  system  can  be  solved  for  different  loads  without  recomputing  the 
coefficients  of  the  system.  This  feature  is  particularly  useful  when  the  constructed  multiple 
crack  solution  is  used  in  the  finite  element  alternating  method,  where  it  is  necessary  to 
evaluate  the  solution  for  the  same  cracks  under  different  loadings  during  the  alternating 
procedure. 
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